//**********************************************************************
//	Programmer	-	Sujal Bista
//	File Name	-	mathFunctions.h
//	Description	:	
//					basic math functions
//
//**********************************************************************

#include "mathFunctions.h"


//computes bezier curve Weights
const float4 bezierWeights(float t)
{
	float temp0=((1-t)*(1-t));
	float temp1=t*t;
	return float4(temp0*(1-t),3*temp0*t,3*(1-t)*temp1,temp1*t);
}

//computes b-spline Weights
const float4 bSplineWeights(float t)
{
	return float4((((-t+3)*t-3)*t+1)/6,(((3*t-6)*t)*t+4)/6,(((-3*t+3)*t+3)*t+1)/6,(t*t*t)/6);
}

//computes Catmull Rom-spline Weights
const float4 catmullRomSplineWeights(float t)
{
	return float4(((-t+2)*t-1)*t/2,(((3*t-5)*t)*t+2)/2,((-3*t+4)*t+1)*t/2,((t-1)*t*t)/2);
}

//computes gaussian Weight
float computeGaussianWeight(const float2 &position,const float2 &center,const float &standardDeviationX,const float &standardDeviationY)
{
	return (1/(2*PI*standardDeviationX*standardDeviationY))*exp(-(((position.x-center.x)*(position.x-center.x))/(2*standardDeviationX*standardDeviationX)+((position.y-center.y)*(position.y-center.y))/(2*standardDeviationY*standardDeviationY)));
}

//computes gaussian kernel
float *computeGaussianKernel(unsigned int radius, float standardDeviationX, float standardDeviationY)
{
	unsigned int size=radius*2;
	float *kernel=new float[size*size];
	float2 center((float)radius,(float)radius);
	
	//compute the gaussian kernel
	for(unsigned int y = 0; y < size; ++y )
	{
		for(unsigned int x = 0; x < size; ++x )
		{
			float2 position((float)y+.5f,(float)x+.5f);
			kernel[y*size+x]=computeGaussianWeight(position,center,standardDeviationX,standardDeviationY);
		}
	}
	return kernel;
}

/////////////////////////////////////////////////////
//functions that effect float4x4S
/////////////////////////////////////////////////////

//computes area of a triangle
float area(const float3 &a,const float3 &b,const float3 &c)
{	
	float3 base=b-a;//base
	float baseLength=length(base);//length of base
	if(baseLength==0)
		return 0;
	float3 ac=c-a;//another edge

	float3 cHat=a+(base/baseLength)*(dot(ac,base)/baseLength); //projection of c on ab
	float3 height=c-cHat;

	return baseLength*length(height)*.5f;
}


//distance from point to line
float distance(const float3 &p0,const float3 &p1,const float3 &pos)
{
	float c1,c2;

	float3 v=p1-p0;
	float3 w=pos-p0;

	//calculate the projectio of p1 onto  the line segment
	c1=dot(w,v);
	if(c1<=0)
		return length(pos-p0);

	c2=dot(v,v);
	if(c2<=c1)
		return length(pos-p1);

	//how far is the projection from c1 [0-1]
	c1=c1/c2;

	//simple paramatic eqn
	w=p0+v*c1;

	return length(pos-w);
}

//distance from point to line
float distance(const float2 &p0,const float2 &p1,const float2 &pos)
{
	float c1,c2;

	float2 v=p1-p0;
	float2 w=pos-p0;

	//calculate the projectio of p1 onto  the line segment
	c1=dot(w,v);
	if(c1<=0)
		return length(pos-p0);

	c2=dot(v,v);
	if(c2<=c1)
		return length(pos-p1);

	//how far is the projection from c1 [0-1]
	c1=c1/c2;

	//simple paramatic eqn
	w=p0+v*c1;

	return length(pos-w);
}

